home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / rng / cmrg.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  5.3 KB  |  198 lines
  1. /* rng/cmrg.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_rng.h>
  23.  
  24. /* This is a combined multiple recursive generator. The sequence is,
  25.  
  26.    z_n = (x_n - y_n) mod m1
  27.  
  28.    where the two underlying generators x and y are,
  29.  
  30.    x_n = (a_{1} x_{n-1} + a_{2} x_{n-2} + a_{3} x_{n-3}) mod m1
  31.    y_n = (b_{1} y_{n-1} + b_{2} y_{n-2} + b_{3} y_{n-3}) mod m2
  32.  
  33.    with coefficients a11 ... a23,
  34.  
  35.    a_{1} = 0,     a_{2} = 63308, a_{3} = -183326
  36.    b_{1} = 86098, b_{2} = 0,     b_{3} = -539608
  37.  
  38.    and moduli m1, m2,
  39.  
  40.    m1 = 2^31 - 1 = 2147483647
  41.    m2 = 2^31 - 2000169 = 2145483479
  42.  
  43.    We initialize the generator with 
  44.  
  45.    x_1 = s_1 MOD m1, x_2 = s_2 MOD m1, x_3 = s_3 MOD m1
  46.    y_1 = s_4 MOD m2, y_2 = s_5 MOD m2, y_3 = s_6 MOD m2
  47.  
  48.    where s_n = (69069 * s_{n-1}) mod 2^32 and s_0 = s is the
  49.    user-supplied seed.
  50.  
  51.    NOTE: According to the paper the initial values for x_n must lie in
  52.    the range 0 <= x_n <= (m1 - 1) and the initial values for y_n must
  53.    lie in the range 0 <= y_n <= (m2 - 1), with at least one non-zero
  54.    value -- our seeding procedure satisfies these constraints.
  55.  
  56.    We then use 7 iterations of the generator to "warm up" the internal
  57.    state.
  58.  
  59.    The theoretical value of z_{10008} is 719452880. The subscript 10008
  60.    means (1) seed the generator with s=1, (2) do the seven warm-up
  61.    iterations that are part of the seeding process, (3) then do 10000
  62.    actual iterations.
  63.  
  64.    The period of this generator is about 2^205.
  65.  
  66.    From: P. L'Ecuyer, "Combined Multiple Recursive Random Number
  67.    Generators," Operations Research, 44, 5 (1996), 816--822.
  68.  
  69.    This is available on the net from L'Ecuyer's home page,
  70.  
  71.    http://www.iro.umontreal.ca/~lecuyer/myftp/papers/combmrg.ps
  72.    ftp://ftp.iro.umontreal.ca/pub/simulation/lecuyer/papers/combmrg.ps */
  73.  
  74. static inline unsigned long int cmrg_get (void *vstate);
  75. static double cmrg_get_double (void *vstate);
  76. static void cmrg_set (void *state, unsigned long int s);
  77.  
  78. static const long int m1 = 2147483647, m2 = 2145483479;
  79.  
  80. static const long int a2 = 63308, qa2 = 33921, ra2 = 12979;
  81. static const long int a3 = -183326, qa3 = 11714, ra3 = 2883;
  82. static const long int b1 = 86098, qb1 = 24919, rb1 = 7417;
  83. static const long int b3 = -539608, qb3 = 3976, rb3 = 2071;
  84.  
  85. typedef struct
  86.   {
  87.     long int x1, x2, x3;    /* first component */
  88.     long int y1, y2, y3;    /* second component */
  89.   }
  90. cmrg_state_t;
  91.  
  92. static inline unsigned long int
  93. cmrg_get (void *vstate)
  94. {
  95.   cmrg_state_t *state = (cmrg_state_t *) vstate;
  96.  
  97.   /* Component 1 */
  98.  
  99.   {
  100.     long int h3 = state->x3 / qa3;
  101.     long int p3 = -a3 * (state->x3 - h3 * qa3) - h3 * ra3;
  102.  
  103.     long int h2 = state->x2 / qa2;
  104.     long int p2 = a2 * (state->x2 - h2 * qa2) - h2 * ra2;
  105.  
  106.     if (p3 < 0)
  107.       p3 += m1;
  108.     if (p2 < 0)
  109.       p2 += m1;
  110.  
  111.     state->x3 = state->x2;
  112.     state->x2 = state->x1;
  113.     state->x1 = p2 - p3;
  114.     if (state->x1 < 0)
  115.       state->x1 += m1;
  116.   }
  117.  
  118.   /* Component 2 */
  119.  
  120.   {
  121.     long int h3 = state->y3 / qb3;
  122.     long int p3 = -b3 * (state->y3 - h3 * qb3) - h3 * rb3;
  123.  
  124.     long int h1 = state->y1 / qb1;
  125.     long int p1 = b1 * (state->y1 - h1 * qb1) - h1 * rb1;
  126.  
  127.     if (p3 < 0)
  128.       p3 += m2;
  129.     if (p1 < 0)
  130.       p1 += m2;
  131.  
  132.     state->y3 = state->y2;
  133.     state->y2 = state->y1;
  134.     state->y1 = p1 - p3;
  135.     if (state->y1 < 0)
  136.       state->y1 += m2;
  137.   }
  138.   
  139.   if (state->x1 < state->y1)
  140.     return (state->x1 - state->y1 + m1);
  141.   else
  142.     return (state->x1 - state->y1);
  143. }
  144.  
  145. static double 
  146. cmrg_get_double (void *vstate)
  147. {
  148.   return cmrg_get (vstate) / 2147483647.0 ;
  149. }
  150.  
  151.  
  152. static void
  153. cmrg_set (void *vstate, unsigned long int s)
  154. {
  155.   /* An entirely adhoc way of seeding! This does **not** come from
  156.      L'Ecuyer et al */
  157.  
  158.   cmrg_state_t *state = (cmrg_state_t *) vstate;
  159.  
  160.   if (s == 0)
  161.     s = 1;    /* default seed is 1 */
  162.  
  163. #define LCG(n) ((69069 * n) & 0xffffffffUL)
  164.   s = LCG (s);
  165.   state->x1 = s % m1;
  166.   s = LCG (s);
  167.   state->x2 = s % m1;
  168.   s = LCG (s);
  169.   state->x3 = s % m1;
  170.  
  171.   s = LCG (s);
  172.   state->y1 = s % m2;
  173.   s = LCG (s);
  174.   state->y2 = s % m2;
  175.   s = LCG (s);
  176.   state->y3 = s % m2;
  177.  
  178.   /* "warm it up" */
  179.   cmrg_get (state);
  180.   cmrg_get (state);
  181.   cmrg_get (state);
  182.   cmrg_get (state);
  183.   cmrg_get (state);
  184.   cmrg_get (state);
  185.   cmrg_get (state);
  186. }
  187.  
  188. static const gsl_rng_type cmrg_type =
  189. {"cmrg",            /* name */
  190.  2147483646,            /* RAND_MAX */
  191.  0,                    /* RAND_MIN */
  192.  sizeof (cmrg_state_t),
  193.  &cmrg_set,
  194.  &cmrg_get,
  195.  &cmrg_get_double};
  196.  
  197. const gsl_rng_type *gsl_rng_cmrg = &cmrg_type;
  198.